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Abstract. This paper summarises an investigation of chaos in a toy potential which 
mimics much of the behaviour observed for the more realistic triaxial generalisations 
of the Dehnen potentials, which have been used to model cuspy triaxial galaxies 
both with and without a supermassive black hole. The potential is the sum of an 
anisotropic harmonic oscillator potential, Vo ~ ^{a'^x'^ + b'^y^ + z'^), and a spherical 
Plummer potential, Vp — —MsH/Vr^ + e^, with = x'^ + y^ + z^. Attention focuses 
on three issues related to the properties of ensembles of chaotic orbits which impact 
on chaotic mixing and the possibility of constructing self-consistent equilibria: (1) 
What fraction of the orbits are chaotic? (2) How sensitive are the chaotic orbits, i.e., 
how large are their largest (short time) Lyapunov exponents? (3) To what extent 
is the motion of chaotic orbits impeded by Arnold webs, i.e., how 'sticky' are the 
chaotic orbits? These questions are explored as functions of the axis ratio a : b : c, 
black hole mass Mbh, softening length e, and energy E with the aims of understand- 
ing how the manifestations of chaos depend on the shape of the system and why 
the black hole generates chaos. The simplicity of the model makes it amenable to a 
perturbative analysis. That it mimics the behaviour of more complicated potentials 
suggests that much of this behaviour should be generic. 
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1. Motivation 

Over the past decade it has become apparent that many early-type 
galaxies are genuinely three-dimensional objects, neither spherical nor 
axisymmetric (c/. Kormendy and Bender 1996) that they typically have 
a central density cusp (c/. Lauer et al 1995) and that the center of the 
galaxy often contains a supermassive black hole (c/. Kormendy and 
Richstone 1995). However, there is substantial numerical evidence (c/. 
Pridman and Merritt 1997) that the combination of triaxiality and a 
cusp typically yields potentials that admit a large amount of chaos; and 
there is also evidence indicating that the introduction of a supermassive 
black hole into a cuspy, triaxial potential causes a further increase 
in the abundance of chaos (c/. Valluri and Merritt 1998). Numerical 
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explorations indicate further (c/. Siopis and Kandrup 2000) that the 
phase space associated with a cuspy triaxial potential, cither with or 
witliout a black hole, can be very complex, laced with an intricate 
Arnold (1964) web that can seriously impede phase space transport, so 
that chaotic orbits in these potentials can be very 'sticky' in the sense 
recognised originally by Contopoulos (1971). 

However, there are at least two important questions which have not 
yet been addressed adequately. Why does the combination of triaxiality 
and a cusp and/or black hole lead to chaos? And to what extent is the 
observed behaviour generic? Almost all work hitherto has focused on 
one class of model potentials, namely the triaxial generalisations of 
the spherical Dehnen (1993) potentials first considered by Merritt and 
Pridman (1996). 

Merritt and Valluri (1996) have argued that a black hole or cusp 
tends to induce chaos because it serves as a near-singular perturbation 
that destabilises what would be regular box orbits if the galaxy had a 
smooth core and there were no black hole. This seems quite reasonable. 
However, Udry and Pfenniger (1988) have noted a fact well known from 
nonlinear dynamics (c/. Lichtenberg and Lieberman 1992), namely that 
chaos often arises simply by breaking a symmetry or combining incom- 
patible symmetries. One might therefore argue that a supermassive 
black hole acts as a source of chaos because of symmetry-breaking: 
the unperturbed galaxy is triaxial, whereas the perturbation manifests 
spherical symmetry. This could, e.g., account for the fact that, as a 
source of chaos, even a comparatively low mass black hole can be more 
important than a steep cusp. As Valluri and Merritt (1998) noted, 'even 
a modest black hole, with a mass ~ 0.3% the mass of the galaxy, is as 
effective as the steepest central density cusp at inducing stochastic 
diffusion.' In any event, one knows that a black hole or cusp does 
not always result in chaos. By carefully maintaining symmetries, one 
can, e.g., construct intcgrablc generalisations of the Stacckcl potentials 
which incorporate a supermassive black hole (Sridhar and Tousma, 
1996; Sridhar and Tousma, 1997). 

One objective here is to address the question of why cusps and black 
holes induce chaos by studying orbits in an extremely simple class of 
potentials which appear to exhibit behaviour qualitatively similar to 
what is observed the triaxial Dehnen potentials with or without a black 
hole. The potential, 



with = + -\- , is the sum of an anisotropic harmonic potential 
and a spherical Plummer potential. There are five parameters which 



y(x,y,z) = -(aV + 5V + cV) 



Mbh 



(1) 



_j_ g2 
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can be varied independently, namely a, b, and c, which determine the 
frequencies and the shape of the unperturbed oscillator, AIbh, the 
mass of the black hole, and the softening length e. That much of the 
behaviour exhibited by orbits in the triaxial Dehnen potentials can be 
reproduced by such a simple model suggests that it may be generic. 
That the model is so simple implies that it is amenable to an analytic 
treatment with the black hole viewed as a perturbation of the integrable 
oscillator potential. 

The numerical component of the analysis focuses on three impor- 
tant characterisations of the chaos: (1) What fraction of the orbits are 
chaotic? (2) How large is the largest Lyapunov exponent associated 
with the chaotic orbits? (3) How 'sticky' are the chaotic orbits? Can 
they travel relatively unimpeded throughout the accessible phase space, 
or is their motion obstructed significantly by the Arnold web? 

The relative abundance of chaotic orbits is important given the 
general presumption that one requires a substantial number of regular, 
or nearly regular, orbits to serve as a skeleton for the equilibrium (c/. 
Binney 1978). The degree of sensitivity exhibited by orbits is important 
because it impacts on the efficacy of chaotic mixing ( cf. Kandrup and 
Mahon 1994, Merritt and Valluri 1999). The existence of "sticky" orbits 
is of potential importance because these orbits can in principal be used 
to support regular structures in phase space regions where few, if any, 
regular orbits exist {cf. Wozniak 1993, Kaufmann and Contopoulos 
1996, Patsis, Athanassoula, and Quillen 1997). 

Determining the answers to these questions as functions of the axis 
ratio a -.b : c and the black hole mass Mbh provides information about 
how the manifestations of chaos vary with the degree of triaxiality and 
the mass of the black hole. Determining the effects of varying e can 
shed insights into why the chaos arises. For sufficiently small softening, 
the exact value of e is immaterial, but when e becomes too large its 
value begins to have an important effect. 

If the only function of the black hole were to act as a near-singular 
perturbation, one might expect that, for ensembles of orbits with fixed 
energy E, unperturbed frequencies a, b, and c, and softening e, the 
relative abundance of chaotic orbits should exhibit a very simple de- 
pendence on Mbh- For Mbh = and Mbh — > oo the potential is 
integrable. However, for other values V is nonintegrable and, for in- 
termediate values far from Mbh = and oo, one might expect large 
measures of chaotic orbits. In particular, one might expect that a plot 
of n{MBH), the relative measure of chaotic orbits, as a function of 
Mbh would exhibit a single maximum at an intermediate value where 
Mbh is neither too small nor too large. This is not what the numerical 
simulations reveal. Rather, one finds that the abundance of chaos varies 
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in a more irregular fashion which, as discussed in Section 2, would 
suggest that chaos arises from a resonance overlap between the natural 
frequencies of the oscillator and the black hole. 

Section 2 of this paper focuses on the statistical properties of short 
time Lyapunov exponents. Section 3 summarises an analysis that fo- 
cused primarily on the Fourier spectra of both individual orbits and 
orbit ensembles. Section 4 concludes by summarising the results of 
Sections 2 and 3 and commenting on their physical implications. 



2. Lyapunov Exponents for Orbit Ensembles 
2.1. WHAT WAS COMPUTED 

Most of the experiments focused on the statistical properties of ensem- 
bles of orbits evolved for comparatively short times. However, these 
were supplemented by very long time integrations of individual initial 
conditions used to extract estimates of the values of the largest Lya- 
punov exponents for fixed values of frequencies a, b, c, mass Mbh, 
softening e, and energy E. Three different questions were addressed: 

1. For fixed a, b, c, MbHi and E, how do things depend on e? For 
sufficiently small e, the precise value is immaterial but for larger e this 
is no longer true. 

2. For fixed a, b, c, e, and E, how do things depend on Mbh'^ The 
most naive picture might suggest that n{MBH), the relative measure 
of chaotic orbits, should vary smoothly with Mbh and that a plot of n 
as a function of Mbh should exhibit no structure. 

3. For fixed e, Mbh, and E, how do things depend on the axis ratio a : 
b : c? Conventional wisdom would suggest that larger deviations from 
spherical symmetry imply more chaos, so that one might expect the 
relative abundance of chaos to increase monotonically with increasing 
asphericity. 

The experiments were performed for frequencies a, b, and c ~ 1 and 
energies £^ ~ 0.1 — 1, this corresponding physically to a galaxy of mass 
M ~ 1, linear size r 1, and characteristical orbital, or dynamical, 
time scale to ~ 2 — 3. The black hole mass was taken in the range 
10-^ <Mbh<1- Smaller black hole masses are not large enough to 
induce significant stochasticity; much larger masses seem completely 
unrealistic. 

For parameters in these general ranges, the exact value of e seems 
unimportant, provided only that < 10~^ or so. For this reason, most 
of the experiments were performed assuming = 10^^, the same value 
used in Siopis &: Kandrup (2000). A 'typical' choice of frequencies was 
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= 1.25, = 1.0, and = 0.75, corresponding to a moderately 
triaxial configuration. 

For fixed frequencies a, b and c, mass Mbh, and energy E, the value 
of the largest Lyapunov exponent was estimated by selecting six or 
more initial conditions corresponding to chaotic orbits and integrating 
each initial condition for a total time t = 256000 while simultane- 
ously tracking the evolution of a small initial perturbation periodically 
renormalised in the usual way (c/. Lichtenberg and Lieberman 1992). 
The relative abundance of chaotic orbits was estimated by sampling the 
phase space to generate a collection of 2000 initial conditions and evolv- 
ing each of these initial conditions, along with a small perturbation, for 
a time t = 2048 which, for the energies considered here, corresponded 
to ~ SOOti). The value of the largest short time Lyapunov exponent 
X {cf. Grassbcrgcr et al 1988) was then used as a diagnostic to deter- 
mine whether or not the initial condition was chaotic. (Integrations for 
shorter times would seem better motivated astronomically, but it was 
found that, for much smaller it was often impossible to determine 
with any reliability whether or not any given segment corresponded 
to a chaotic orbit. Indeed, even for integrations for times as long as 
t = 2048 it is possible that some very sticky chaotic segments were 
misinterpreted as regular segments!) The initial conditions were gener- 
ated as in Siopis and Kandrup (2000), using a scheme which generalises 
Schwarzschild's (1993) algorithm to construct a representative library 
of orbits. 

The resulting short time x's were examined to determine three 
things, namely (1) the relative number of chaotic orbits; (2) the de- 
gree to which different chaotic orbit segments with the same energy 
do, or do not, have the same short time x's, and (3) the degree to 
which all the chaotic segments are characterised by short time x's well 
separated from zero. The behaviour of orbit segments over shorter, 
more astronomically relevant, time scales was addressed by partitioning 
each of the 2000 segments of length t = 2048 into eight segments of 
length t = 256 and extracting a short time Lyapunov exponent for each 
segment {cf. Kandrup and Mahon 1994). The resulting 16000 segments 
were analysed to determine a distribution of short time Lyapunov ex- 
ponents, N[x{t = 256)], which probes the varying degrees of chaos 
exhibited by a representative collection of orbits, both regular and 
chaotic, over an interval of time comparable to the age of the Universe. 
Finally, as another probe of the total amount of chaos, independent of 
distinctions amongst regular, sticky, and wildly chaotic orbits, the out- 
put was analysed to extract the average short time Lyapunov exponent 
for all the orbits in each ensemble. 
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Most of the experiments focused on three representative energies, 
E = 1.0, 0.6, and 0.4. A variety of axis ratios were considered. One set 
of experiments involved considering frequencies a, b, and c for which 

: 5^ : = 1 + A : 1 : 1 — A and exploring the effects of increasing 

A. 

2.2. WHAT WAS FOUND 

The black hole acts as an important source of chaos whenever con- 
ditions are such that the force that it exerts on a test particle often 
becomes comparable to the force associated with the oscillator. For 
small e, a mass Mbh > ^min ~ 10~^ appears required to induce an 
appreciable amount of obvious chaos. The precise value of e seems 
unimportant provided only that e is very small compared with the char- 
acteristic size of a stellar orbit. If the black hole mass be frozen at some 
value Mbh significantly above Mmin and e increased systematically, 
chaos remains common until the maximum black hole force ~ Mbh / ^ 
becomes smaller than a characteristic value somewhat less than unity, 
this corresponding to the typical force exerted by the oscillator. It thus 
appears that it is not the singular character of the perturbation per se 
that is responsible for triggering chaos; rather what is important is that 
the force associated with the black hole can become comparable to, or 
larger than, the force associated with the unperturbed oscillator. In 
other words, mixing two integrable potentials of comparable strength 
is responsible for the chaos. 

One example of how the amount of chaos depends on e is provided 
by Figure 1, each panel of which exhibits the computed short time 
x(t = 2048) for a collection of 2000 initial conditions evolved in the 
potential (1) with c? = 1.25, \? = 1.0, (? = 0.75, and Mbh = 0.1 
with fixed energy E = 1.0. Each panel was constructed by computing 
2000 short time x's, ordering the orbits in terms of the computed x, 
and then (c/. Siopis and Kandrup 2000) plotting the resulting x's as a 
function of particle number from 1 to 2000. The six panels exhibit data 
for values of e ranging from = 10~^ to 10~^. The computed x's for 

= 10~^, 10~^, and 10~^ are virtually identical. For = 10~^ chaos is 
weaker in two senses, namely that the relative number of chaotic orbits 
is smaller and that the typical x associated with a chaotic orbit has 
also decreased. For = 10~^, chaos is almost completely suppressed. 
For = 10~^, there are a few orbits with somewhat larger x's than 
was observed for = 10~^, 10~^, and 10~^, each corresponding to a 
particle which passed extremely close to the black hole. 

Overall, chaotic orbits tend to be extremely sticky. In particular, 
even for t = SOOto or larger, the computed values of x need not man- 
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Figure 1. Short time = 2048) for 2000 orbits evolved in the potential (1) with 
a? = 1.25, &^ = 1.0, = 0.75, E = 1.0, and Mbh = 0.1 with variable e. (a) 
e = 10-'^ ^ (b) e = 10-1^ (c) e = IQ-^ ^ (d) e = lO-' ^. (e) e = IQ-^'^ . (f) 
e = 10-3°. 



ifest a clear and unambiguous separation between regular and chaotic 
behaviour. If, as in Figure 1, an ordered list of x's is plotted as a 
function of particle number, there is in general no clean break between 
the regular and chaotic orbits. Oftentimes the only way to infer a 
reasonable demarcation between regular and chaotic behaviour is to 
look for a 'kink' in the slope. However, the relative abundance of regu- 
lar/nearly regular as opposed to wildly chaotic orbits can be extracted 
by examining distributions of short time Lyapunov exponents, N[x{f^]- 
Typical examples of such distributions are provided in Figures 2 and 3, 
which exhibit N[x{t)] for t = 256 and t = 2048 for the 2000 orbit 
ensembles evolved with = 1.25, 6^ = 1.0, (? = 0.75, e = 10~^, and 
E = 1.0 with variable black hole masses ranging between Mbh = 10~^-^ 
and Mbh = 1- 

In each of these panels, the data were so plotted that the horizontal 
axis terminates at the right at a value Xmax corresponding to the 
largest value of x{t = 256) assumed by any of the 16000 segments. 
It is thus evident that, for small Mbh, there are a very few orbits 
with comparatively large values of x, even though the overwhelming 
majority of the orbits have much smaller values. For relatively large 
values of Mbh, e.g., Mbh = 10~°'^^, there are comparable numbers 
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Figure 2. Distribution of short time Lyapunov exponents, N[x(t = 256)] (thick 
hnes) and N[xit = 2048)] (thin hnes) generated from 2000 orbit segments evolved 
in the potential (1) with = 1.0, 6^ = 1.25, = 0.75, E = 1.0, e = 10"^, and 
variable Mbh- (a) logjo Mbh = -3.5. (b) log^g Mbh = -3.0. (c) log^j Mbh = -2.5 
.(d) logio Mbh = -2.25. (e) logjo Mbh = -2.0. (f) logio Mbh = -1.75. 



of both regular /nearly regular and wildly chaotic orbits. For somewhat 
lower values, the abundance of wildly chaotic orbits tends to decrease 
with increasing Mbh, although the trend is not completely uniform. 
Thus, e.g., there are more wildly chaotic orbits for Mbh = 10"'^'^ and 
Mbh = 10-^-^ than for Mbh = 10"^ 

For small values of Mbh, it becomes exceedingly difficult to estimate 
the fraction of the orbits which are truly chaotic. In particular, for 
Mbh = 10~^ and less, the distribution N[x\ manifests a single sharp 
peak at a value close to zero. This might suggest that the orbits are all 
regular; but N[x] also has a tail extending to much larger values which 
is clearly indicative of chaos. One possible interpretation is that, for 
small Mbh, many/all of the orbits are in fact chaotic, but that most of 
these chaotic orbits behave in a near-regular fashion almost all of the 
time. For Mbh = 0, each initial condition corresponds to a regular box 
orbit which densely fills a region in configuration space that includes 
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Figure 3. Distribution of short time Lyapunov exponents, N[x(t = 256)] (thick 
hnes) and N[x{t = 2048)] (thin hnes) generated from 2000 orbit segments evolved in 
the potential (1) with = 1.0, &^ = 1.25, = 0.75, E = 1.0, e = 10"^, and variable 
Mbh- (a) log^o Mbh = -1.5. (b) log^g Mbh = -1.0. (c) log^o Mbh = -0.75. (d) 
logio Mbh = -0.5. (e) logm Mbh = -0.25. (f) log^ Mbh = 0. 



the origin. One might perhaps expect that, for sufficiently small Mbh, 
the orbits remain "nearly boxy" and again follow trajectories which 
occasionally bring them arbitrarily close to the origin. In the limit 
e ^ 0, one would thus anticipate that every orbit eventually experi- 
ences perturbations of arbitrarily large amplitude which trigger chaotic 
behaviour. Making e finite introduces a bound to the perturbing force 
but, for e as small as 10~^, the value used here, this still corresponds 
to a maximum force of amplitude lO'^ Mbh- 

This interpretation is corroborated by the fact that similar distri- 
butions of short time Lyapunov exponents can be generated from a 
single initial condition integrated for very long times. Specifically, if a 
single initial condition with the same energy be integrated for a time 
k X 256 and short time exponents computed separately for each t = 256 
segment, the resulting N[x] is often nearly indistinguishable from the 
N[x] generated from an ensemble of initial conditions. This is, e.g., 
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Figure 4- (a) Distribution of short time Lyapunov exponents, N[x{t — 256)] gener- 
ated from a single initial condition with = 1.0, 6^ = 1.25, = 0.75, e = 10~^, 
E = 1.0, and Mbh = integrated for a total time t = 8000 x 256. The initial 
condition corresponded to an orbit which for early times was wildly chaotic, (b) 
N[x{t = 256)] for another initial condition with the same parameter values, in this 
case corresponding at early times to a 'sticky,' near-regular orbit. 



evident from Figure 4, which exhibits distributions generated for two 
different initial conditions with E = 1.0 integrated with = 1.25, 
6^ = 1.0, = 0.75, e = 10"^, and Mbh = 10"^. Each initial condition 
was integrated for a total time t = 8000 x 256 with energy conserved 
to better than one part in 10^. One corresponded at early times to a 
wildly chaotic orbit with = 256) ~ 0.06; the other corresponded 
initially to a very 'sticky,' near-regular orbit. 

Figure 5 exhibits estimates of x, the largest Lyapunov exponent, and 
n{MBH), the relative fraction of obviously chaotic orbits, as functions of 
Mbh for representative ensembles computed with = 1.25, 6^ = 1.0, 

= 0.75, and e = 10~^ for three different energies, namely E = 1.0, 
E = 0.6, and E = 0.4. Only the data for Mbh > 10"^'^ are exhibited 
since, as noted above, for smaller values it becomes extremely difficult 
to determine whether any given orbit is regular or chaotic. However, 
it does seem clear that there is little or no chaos when the black hole 
mass is as small as Mbh = 10^^, and that there is a comparatively 
rapid increase in the importance of chaos for Mbh ~ 10~^. This rapid 
increase is probably not generic. Rather, as discussed in Section 4, the 
fact that chaos appears to 'turn on' abruptly likely reflects the fact 
that, for this potential, all the unperturbed {Mbh = 0) orbits with 
fixed energy E have exactly the same natural frequencies and densely 
sample the same region in configuration space, so that 'if you've seen 
one orbit you've seen them all.' 

If the black hole mass be increased from Mbh = 10^^'^ to larger 
values, the overall abundance of chaos tends to decrease. However, this 
decrease is not completely uniform. For all three energies, one observes 
a pronounced spike in the distribution n{MBH) where the relative 
measure of chaotic orbits increases abruptly. The existence of these 
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Figure 5. (a) Estimates of the largest Lyapuuov exponent for chaotic orbits evolved 
in the potential with = 1.0, 6^ = 1.25, = 0.75, e = 10"^, and E = 1.0 as a 
function of Mbh- (b) The fraction of chaotic orbits from a representative ensemble 
with the same a, 6, c, e, and £ as a function of Mbh- (c) and (d) The same as (a) 
and (b) but for E = 0.6. (e) and (f) The same as (a) and (b) for E = 0.4. 



spikes indicates that the relative abundanee of ehaos does not exhibit 
a completely trivial dependence on Mbh- The fact that the locations 
of these spikes varies smoothly with energy suggests strongly that this 
nontrivial variation is a resonance phenomenon. 

Similarly, for fixed a, b, c, e, and E, the largest Lyapunov exponent, 
as defined in an asymptotic t — > oo limit, does not exhibit a completely 
trivial dependence on Mbh- In particular, x is not a monotonically 
increasing function of Mbh- However, for ensembles integrated for 
times t ~ 100 — SOOiu, it appears that the largest short time exponents 
observed for any of the orbits do increase monotonically with increasing 
Mbh- The fact that, for certain ranges of Mbh, the asymptotic % 
decreases reflects the fact that, for these values of Mbh, it is especially 
likely for a chaotic orbit to be 'nearly regular' for long times, during 
which the short time x is comparatively small. That the locations 
of these 'dips' varies smoothly with energy again suggests that some 
resonance phenomenon has come into play. 

Overall, chaos seems to become more pronounced for larger devia- 
tions from axisymmetry but, once again, the trend is not completely 
uniform. Rather, it is again possible to observe dips and/or spikes in the 



CM.tex; 1/02/2008; 20:42; p. 11 



12 



H. E. Kandrup and I. V. Sideris 




Figure 6. (a) The average short time Lyapunov exponents, (x) (diamonds), and the 

associated dispersions, (triangles), for representative ensembles with e = 10"'^, 
E = 1.0, Mbh = 0.1, and : fe^ : = 1 + A : 1 : 1 - A for variable A. (b) The 
same for an ensemble with E = 0.6. (c) The fraction n of wildly chaotic orbits for 
the E = l.Q ensemble in (a), (d) The same for the E = 0.6 ensemble. 

relative measure of chaotic orbits. This is, e.g., ihustrated by Figure 6, 
which exhibits (x), the mean value of the largest Lyapunov exponent, 
C7^, the associated dispersion, and n, the relative abundance of strongly 
chaotic orbits for sequences of models with fixed E, Mbh = 0.1, 
e = 10-2, and : 6^ : = 1 + A : 1 : 1 - A with < A < 0.8. 
The two left panels were computed for ensembles with E = 1.0; the 
two right panels were computed for E = 0.6. 

3. Spectral Properties of Orbits and Orbit Ensembles 
3.1. WHAT WAS COMPUTED 

The experiments described here aimed to provide additional insights 
into the origins of chaos in the potential (1) by examining the Fourier 
spectra of representative orbits and orbit ensembles; testing the validity 
of a perturbative analysis which views Vp = —Mbh I ^ r"^ + ^ as a 
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perturbation of the unperturbed oscillator potential; and examining 
how, for individual orbits, the degree of instability exhibited at different 
points correlates with the location of the orbit in configuration space. 

Two classes of experiments were considered. The first involved com- 
puting Fourier spectra for each orbit in several of the ensembles de- 
scribed in the preceding section. This facilitated an improved under- 
standing of how various orbits with the same E, a, b, c, e, and Mbh 
can differ one from another. 

The other class of experiments explored the question of how, for a 
fixed set of initial conditions, evolved with the same choices of a, b, 
c, and e, varying the central mass Mbh impacts bulk statistical prop- 
erties. This entailed selecting fifty representative phase space points, 
uniformly sampling the interval 0.16 < r < 1.4, with Vy = Vz = 
and Vx so chosen that E = 1.0 for Mbh = 0.1; and then evolving 
these into the future for a time t = 2048 with = 1.25, 6^ = 1.0, 

= 0.75, and e = 10~^, varying the black hole mass incrementally from 
Mbh = to Mbh = 0.4 in steps of 0.01. An analysis of the resulting 
Fourier spectra facilitated a quantitative characterisation of the extent 
to which perturbation theory correctly predicts changes in frequencies. 
A comparison of short time Lyapunov exponents and properties of 
the Fourier spectra for the orbits corroborated the expectation that, 
as for other two- and three-dimensional potentials (c/. Kandrup et al 
1997, Siopis et al 1998), orbit segments with larger Lyapunov exponents 
typically have more 'complex' Fourier spectra, with substantial power 
at a larger number of different frequencies. 

The notion of 'complexity' exploited here differs slightly from that 
adopted in earlier papers. In those papers, the complexity n{k) was 
defined as the number of frequencies required in a discrete Fourier 
spectrum to capture a fixed fraction k of the total power. Here the 
complexity was defined instead as equaling the number of frequencies 
which contained power greater than or equal to a fixed fraction k of 
the power at the peak frequency. The results exhibited here involved 
the choice A; = 0.7 and a complexity 

= ^ [nx{i^) + ny{uj) + nx{uj)] (2) 

where, e.g., Hx denotes the number of frequencies for motion in the x- 
direction. The alternative prescription exploited here reflects the fact 
that there is no obvious unique definition of complexity, and that 
different definitions yield comparable results. 

The analytic computation of modified frequencies was effected us- 
ing standard second order perturbation theory (c/. Lichtenberg and 
Lieberman 1992). Relevant formulae are summarised in an Appendix. 



CM.tex; 1/02/2008; 20:42; p. 13 



14 



H. E. Kandrup and I. V. Sideris 



10 
8 

6 

c 

4 

2 



0.00 0.10 0.20 0.30 0.40 
Mbh 

Figure 7. n{u3), the mean number of frequencies with amplitude greater than 0.7 
times the frequency with peak power (sohd curve), and Xi the mean short time 
Lyapunov exponent (dashed curve), computed for the same set of fifty different 
initial conditions with E = 1.0 evolved in the potential (1) for t — 2048 with 
= 1.25, = 1.0, and = 0.75, for a range of black hole masses < Mbh < 0.4. 

3.2. WHAT WAS FOUND 

Overall, the complexities of different orbit segments correlate well with 
the values of the largest short time Lyapunov exponent x- orbits for 
which power is spread over a larger number of frequencies tend sys- 
tematically to have larger values of x- That this should be the case for 
different orbits evolved in the same potential with the same energy is 
hardly surprising, given results for other three-dimensional potentials 
(c/. Siopis et al 1998). Less obvious, however, but also true is that such 
a correlation persists even when comparing ensembles with different 
values of Mbh- This is, e.g., evident from Figure 7, which plots as a 
function of Mbh the average values of n{uj) (solid curve) and x (dashed 
curve) for the aforementioned fifty orbit ensembles. It is clear that the 
two plots are quite similar, each peaking at a mass Mbh ~ 0.2 where, 
consistent with Figures 2 and 3, the relative measure of wildly chaotic 
orbits is especially high. 

For most values of Mbh there are relatively few orbits with very 
large values of x- uiost of the orbits are either regular or nearly reg- 
ular. For example, for E = 1.0, = 1.25, 6^ = 1.0, and = 0.75, 
there arc large measures of wildly chaotic orbits only for M bh between 
about 10~^'^ and lO"'''^. For much larger or smaller values Mbh one 
finds that N[x\, the distribution of short time Lyapunov exponents, 
is a decreasing function of x- Thus, it is hardly surprising that A^[n], 
the distribution of complexities, should be a decreasing function of n. 
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Figure 8. The relative number of orbits with different complexities n(ui), computed 
for the same ensemble used to generate Figure 3 b. 



Less obvious, however, but also true is that even for a coUection of 
wildly chaotic orbits, where N[x\ peaks at some intermediate value 
of X: N[n] can be a sharply decreasing function of n. This is, e.g.^ 
evident from Figure 8, which plots N\n\ for the same orbit ensemble 
used to generate the distribution N[x\ in Figure 3 b. In this plot, the 
orbits with n{uj) < 2 are almost all regular or nearly regular, whereas 
the orbits with n{uj) > 2 are almost all wildly chaotic. The obvious 
point is that, for n{Lo) > 2, the distribution N[n] decreases in a near- 
uniform fashion, even though the corresponding distribution N[x] has 
a pronounced peak at x ~ 0.06. 

Second order perturbation theory can be expected to work reason- 
ably well only when Mbh is relatively small; and, for this reason, 
comparisons between analytic predictions and numerical computations 
were only effected for Mbh = 0.4 and less. Even for such small val- 
ues, however, subtleties remain. By construction, perturbation theory 
computes shifts in three basic frequencies which, for Mbh = 0, reduce 
to a; = a, b, and c. However, one discovers oftentimes that, even for 
regular orbits, Mbh > implies spectra with power at more than three 
frequencies; and, for fully chaotic orbits, the spectra should, strictly 
speaking, be continuous (c/. Tabor 1989), although most of the power 
may be concentrated at or near a small number of frequencies. For this 
reason, identifying the relative error between the analytics and numer- 
ical computations requires a concrete prescription to identify what one 
really means by the 'principal' frequencies for oscillations in the x, y, 
and z directions. As a practical matter, a principal frequency like ujx 
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Figure 9. The fractional error between the values of the 'principal' frequency uix 
computed numerically and calculated analytically for a single orbit with E = l.Q 
evolved in the potential (1) for t = 2048 with = 1.25, = 1.0, and = 0.75 for 
a range of black hole masses < Mbh < 0.4. Squares denote chaotic orbits, whereas 
diamonds denote regular or nearly regular orbits. 

was extracted by identifying those frequencies in the Fourier spectrum 
which have at least 70% as much power as that contained in the fre- 
quency with maximum power, and then computing a power-weighted 
average value for these frequencies. A net error between analytics and 
numerics for a single orbit was derived by averaging over the error 
associated with the three different frequencies, Ux, ^^y, and tOz- 

Overall, it was found that, for Mbh < 0.2 or so, the typical dis- 
crepancy between the analytics and numerics is quite small, less than 
or of order 5%. However, as illustrated in Figure 9, which exhibits 
the error for Ux for a single representative orbit, the typical errors 
become substantially larger for Mbh > 0.2 or so. One might suppose 
that analytic predictions of the 'principal frequencies' for chaotic orbits 
would tend to be much worse than for regular orbits, and that the 
abrupt increase in relative error for Mbh > 0.25 reflects the onset of 
chaos. This, however, is not the whole story. For this initial condition, 
almost all the orbits with Mbh > 0.08 are obviously chaotic, but the 
fractional errors are all less than 0.06 for Mbh < 0.3. 

Nevertheless, it is apparent that, for fixed value of Mbh, the typical 
error grows systematically with increasing complexity, which is hardly 
surprising in light of the fact that, for large complexity, the 'principal' 
frequency really involves an average over several different frequencies. 
This is, e.g., evident in the first panel of Figure 10, which exhibits 
the average error as a function of n{uj) for the same orbit ensemble 
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Figure 10. (a) The average fractional error as a function of n{u>), the number of 
frequencies with amplitude greater than 0.7 times the frequency with peak power 
for the orbit ensemble used to generate Figure 7. The error bars reflect the spread 
of values observed for orbits with the same n{uj). (b) The relative fraction of the 
orbits in the ensemble with different fractional errors. 



used to generate Figure 7. That perturbation theory is more successful 
overall in accounting for regular and nearly regular orbit segments 
than for wildly chaotic segments is evident from the second panel in 
Figure 10, which plots N[e\, the distribution of relative errors. That 
the distribution is sharply bimodal is consistent with the fact that, as 
regards the success of perturbation theory, the orbits divide into two 
distinct populations, one, comprised mostly of regular orbits, for which 
perturbation theory is quite reliable, and another, comprised mostly of 
wildly chaotic orbits, for which perturbation theory works somewhat 
less well. 

But precisely why does the observed chaos actually arise? A detailed 
examination of the local stability of individual orbits provides a simple 
and compelling answer. When an orbit is far from the center, the central 
mass Mbh has only a comparatively minimal effect and the orbit is 
nearly regular. However, when the orbit comes too close to the center 
and Mbh contributes significantly to the total potential experienced 
by the orbit, the combination of the incompatible symmetries leads to 
a resonance overlap and a significant exponential instability. This is, 
e.^., illustrated in Figure 11, which exhibits a representative piece of a 
wildly chaotic orbit segment with E = 0.75 evolved in the potential 
(1) with = 1.25, 6^ = i.o, = 0.75, and Mbh = 0.15. Here 
the solid curve in the top panel displays the phase space separation 
between the original orbit and a perturbed orbit displaced 
originally by a distance \SZ\ = 10~^ and renormalised in the usual 
way (c/. Lichtenberg and Lieberman 1992) at intervals St = 12.0. The 
dashed curve shows an analogous plot of \dZ\ for a regular orbit. The 
bottom panel exhibits r{t), the distance from the origin (solid curve), 
and \Vp\ = Mbh/Vt^ + e'^, the magnitude of the central mass per- 
turbation. The obvious point is that most of the time the perturbed 



CM.tex; 1/02/2008; 20:42; p. 17 



18 



H. E. Kandrup and I. V. Sideris 




Figure 11. (a) The solid curve shows the magnitude of the phase space 

distance between a perturbed and unperturbed chaotic orbit evolved in the potential 
(1) with = 1.25, 6^ = 1.0, = 0.75, Mbh = 0.15, and E = 0.75. For comparison, 
the dashed line shows \5Z\ for a regular orbit. The orbit was rcnormalisod at intervals 
5t = 12.0 at points indicated by vertical stripes in the plot, (b) The solid curve shows 
r(t), the distance from the origin, for the same orbit at the same times. The dotted 
line shows the magnitude of the perturbation, = MbhI^/t^ + e^- The vertical 
line corresponds to a value r = 0.17, for which |T4>| ^ 0.88. 



and unperturbed orbit remain very close together, with httle if any 

systematic exponential divergence; but that when r < 0.17 or so, so 
that \Vp\ becomes as large as ~ 0.88, the two orbits tend to diverge 
significantly. 

In this context, it should be noted that, even for the orbits which 
have been characterised as 'wildly chaotic,' the computed Lyapunov 
exponents are comparatively small when expressed in units oi to- For 
'generic' chaotic potentials, one discovers typically that the largest Lya- 
punov exponent typically assumes a value ~ 0.5tjj^ — l.Ot]^^ which, for 
t£) r-u 2 — 3, implies values of x in excess of 0.2 or so. However, only for 
very large central masses, Mbh > 0.3 or so, did any of the computed 
orbit segments have values of x as large as 0.2. Even the so-called wildly 
chaotic orbits tend to be comparatively regular most of the time, the 
observed chaos resulting from those comparatively rare intervals when 
the orbits come comparatively close to the central mass. 



4. Discussion 



The stickiness manifested by chaotic orbits in the potential explored 
in this paper is strikingly reminiscent of what has been observed for 
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the triaxial analogues of the Dehnen potentials (c/. Merritt and Valluri 
1996, Siopis and Kandrup 2000), especially for those orbits which, in the 
absence of the cusp and/or black hole, would correspond to box orbits. 
Introducing a cusp and/or a supermassive black hole into an otherwise 
smooth triaxial potential can generate a considerable amount of chaos, 
but individual orbit segments can be very nearly regular for very long 
times. Stickiness appears to be a more important phenomenon in cuspy 
triaxial potentials, both with and without a central black hole, than in 
many other chaotic three-dimensional potentials. 

The numerical experiments described in Section 2 - especially the 
fact that quantities like the relative measure of chaotic, as opposed to 
regular, orbits can manifest a nontrivial dependence on both Mbh and 
axis ratio - suggests strongly that the chaos associated with this poten- 
tial results from a resonance overlap between the natural frequencies 
of the unperturbed oscillator and the natural frequencies of the central 
black hole. An examination of the local stability of individual orbits 
indicates that, for intermediate black hole masses, the orbits are only 
very unstable when they are sufficiently close to the black hole that | Vg \ 
and \ Vp\ become comparable in magnitude. The computations for this 
simple potential thus corroborate the conclusions of Valluri and Merritt 
(1998) based on a detailed investigation of the spectral properties of 
orbits in the triaxial Dehnen potentials. However, it is also evident 
that the chaos arises in response to a combination of incompatible 
symmetries: the Plummer and oscillator potentials are each integrable 
separately but, when combined with comparable magnitudes, can yield 
strongly chaotic behaviour. 

The simplicity of the potential (1) was important because it made 
possible the analytic calculations described in Section 3. However, this 
potential is artificial in the sense that, for fixed energy E, every orbit 
unperturbed by a central black hole has the same natural frequencies, 
a fact that renders the conclusions of the paper anomalous in two 
significant respects: (1) As is clear, e.g., from Figures 2 and 3, as Mbh 
increases there is an extremely abrupt transition from a phase space 
hypersurface with essentially no chaos to a phase space that is almost 
completely chaotic. If the harmonic oscillator component of the poten- 
tial is replaced by an anharmonic oscillator, e.g., by allowing for quartic 
corrections, the natural frequencies of different unperturbed orbits with 
the same energy are no longer identical, and one would anticipate the 
transition from little chaos to much chaos to become more gradual. 
(2) If the oscillator is made anharmonic, the nontrivial structures ex- 
hibited by u^Mbh) and n(A) might also be expected to disappear. 
Because an ensemble of fixed energy now involves unperturbed orbits 
with different frequencies, a plot of n{MBH) can be understood as 
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Figure 12. (a) The average short time Lyapunov exponents, (x) (diamonds), and the 
associated dispersions, (triangles), for representative ensembles with e = 10~^, 
E = 1.0, and ai'' :}?■.(? = 1.25 : 1.00 : 0.75 evolved in the anisotropic potential (3) 
with Q = 1.0 and variable Mbh- (b) The average short time Lyapunov exponents, 
(x) (diamonds), and the associated dispersions, a-)^ (triangles), for representative 
ensembles with e = 10"^, E = 1.0, Mbh = 0.1, and : 6^ : = 1 + A : 1 : 1 - A 
for variable A. 



involving a nontrivial superposition of orbits with different unperturbed 
frequencies. 

These expectations were confirmed by repeating the computations 
described in this paper for orbits in the generaUsed potential 

F(x,y,^) = i(aV + 6V + cV) 



+^(aV + 6V + cV)-^, (3) 

for a variety of different values of a. Data for the representative value 
a = I are exhibited in Figure 12. Here the first panel exhibits the 
mean (x) and dispersion a-^ for 2000 orbit ensembles with E = 1.0, 
a2 = 1.25, 6^ = 1.0, c2 = 0.75 and = lO"'' for different values of 
Mbh ■ The second panel exhibits the same quantities for ensembles with 
E = 1.0, Mbh = 0.1, = lO""^, and : ^2 • ^2 = l + A : 1 : 1 - A 
for different values of A. It is evident from panel (a) that that the 
transition to chaos is more gradual than was observed in the harmonic 
potential. Moreover, only for black hole masses somewhat smaller than 
Mbh = lO""^ do ah the orbits seem completely regular. And similarly, 
it is clear that (x) and a-^ vary smoothly with Mbh and A, the resonant 
behaviour observed for a = having washed out because the ensembles 
now involve unperturbed orbits with a variety of frequencies. 
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Appendix 



The Hamiltonian appropriate for motion in the potential (1) can be 
written in the form H = Hq + Hi , where 



is integrable and 

Mbh 

rl\ = — 



(4) 



(5) 



may be viewed as a perturbation that breaks integrabihty. Because 
the integrable component Hq is separable, it is easy to compute the 
generating functional W{Qi, Ji,Q2, J2,Q3, J3), which can be written 
as 



(6) 



where Jj denote the actions and the index i ranges over x, y, and z. 
Given W, it is easy to relate the original canonical variables to the 
actions Jj and a corresponding set of angles 9i, in terms of which the 
Hamiltonian reduces to H = Hq + Hi, where 



Ho = Ji + J2 + J3 



(7) 



and 



Hi = -Mbhx 



^ sin2 AiOi + ^ sin2 A^O^ + ^ sin^ ^3^3 + 
A2 A^ 



-1/2 



(8) 



At this stage, one needs to perform a triple integral, averaging each of 
the three angles over the range < ^ < 2it. Unfortunately, however, 
this is difficult to do in closed form. Nevertheless, one can proceed 
perturbatively by expanding Hi in a Taylor series around some point 
{61,62,93) = {qi,Q2,a3). The result of such a computation is that the 
perturbed frequencies Ui satisfy 



A. 



I oJi 



Hi{MBH,qi,q2,q3,Ji,J2,J3) + 
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It is natural to choose the qi's at or near values for which A^qi = ±7r/2, 
for which values the magnitude of Hi is minimised. 
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